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There is a large variety of turbulence models available. These 
models include direct numerical simulation, large eddy simulation, 
Reynolds stress/ flux model, zero equation model, one equation model , 
two equation k — € model, multiple - scale model, etc. Each turbulence 
model contains different physical assumptions and requirements. The 
natures of turbulence are randomness, irregularity, diffusivity and 
dissipation. In this study, the capabilities of the turbulence 
models, including physical strength, weakness, limitations, as well 
as numerical and computational considerations, have been reviewed. 
Recommendations are made for the potential application of a 
turbulence model in thrust chamber and performance prediction 
programs. In this study, the full Reynolds stress model is 
reco mm ended - In a workshop, specifically called for the assessment 
of turbulence models for application in liquid rocket thrust 
chambers, most of the experts present were also in favor of the 
recommendation of the Reynolds stress model. 
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I. Introduction 


Turbulence consists of random velocity fluctuations, so in 
principle, it must be treated with statistical methods. Turbulent 
momentum fluxes and turbulent diffusion fluxes are usually much 
greater than the mean molecular length scale in the turbulent regime. 
However, the turbulent length scale of the excited modes still remains 
within the band of computationally resolved grids (either in the 
physical or spectral domain) . This restriction imposes a severe 
limitation on the development of the direct numerical simulation 
method representing turbulent flows even with the advancement of 
present day supercomputer technology, in addition to the 
consideration of the mathematical description of reactive flow 
transport due to the randomness in the fields of velocity and 
thermochemical fluctuations (Givi, 1989; 1990). In other words, 
assuming that the direct computation of the usual thermochemical 
balance equations can be carried out with the help of much better 
supercomputers, it is still impossible in practice to accomplish the 
computation of every realization of this ensemble of flow fields, even 
if we are satisfied simply with a brute numerical method (Borghi, 
1988) . 

The nature of turbulence covers a wide spectrum of turbulent 
energy transports in which turbulent energy production is associated 
with smaller wave numbers; while turbulent energy dissipation is 
attached to larger wave numbers. In other words, turbulence consists 
of fluctuations of motions varying continuously with a spectrum of 
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length and time scales. It is obvious that any turbulence closure 
model with a single length scale can not properly cover the transition 
of turbulent energy from the production end to the dissipation end. 

Thermochemical and dynamical flow fields of liquid rocket thrust 
chambers are highly irregular and random in nature. Proper selection 
of the best turbulence model is the key to understanding the 
characterisitics of reacting flow fields. 

Turbulence models available include model-free simulations 
(including large eddy simulation) , which is a direct numerical 
simulation of thermochemical flow fields, and turbulence closure 
models which apply Favre and Reynolds averages on flow fields with 
random fluctuations. Turbulence closure models can be further 
divided into the models which do and do not adopt the concept of 
Boussinesq / s expression. By adopting the concept of Boussinesg's 
expression, Reynold stress and diffusion fluxes become eddy viscosity 
and eddy diffusion fluxes repectively, which assume that the flow 
fields are in the form of isotropic turbulence. This approach has 
produced formulations which are known as: zero equation, one 
equation, two equation (including k— € model) , and multi— scale 
turbulence models. Reynolds stress/flux models were developed by 
the direct derivation of Reynolds stress/diffusion fluxes from the 
mean value balance equations of continuity, momentum, energy and 
species concentration, instead of using eddy viscosity/eddy 
diffusion flux concepts from Boussinesq's expression. Reynolds 
stress/flux models also include algebraic stress/flux and multi- 
scale turbulence. 
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With this large variety of turbulence models available, the 
present report discusses the advantages and disadvantages of applying 
these models in the reacting flow calculation, in particular the 
suitability for random flow field computations within a liquid rocket 
thrust chamber. 

Please note that some text in this report is repetitive for the 
sole purpose of convenience to the reader and to emphasize the 
importance of the topic in question. 


II. Formulations of Turbulent FlPWS In Exacting Fluid s 


It is now common to consider a turbulent flow as an ensemble of 
random flow fields. This includes the following considerations: 
(1) each of these fields does satisfy the classical thermodynamical 
balance equations with a particular set of initial and boundary 
conditions; (2) each field is different due to the differences in 
initial and boundary conditions; and (3) each flow can be very 
different because the sensitivity to boundary conditions is quite 
large. This Turbulent Regime occurs very often when a characteristic 
Reynolds number (a ratio of convective force to viscous force) is 
large enough. Other types of instabilities leading to Tutbulencg 
which are characterized by other numbers could occur, especially for 
reacting flows where highly non-linear terms appear in the equations. 

Assuming that the thermochemical balance equations are valid, it 
is still impossible to compute these flow fields even if we settle for 
a brute numerical method (Borghi, 1988) . The time scales and length 
scales, that we know to exist within the turbulent regime of reacting 
flows, are so small with respect to the time or length scales of 
interest that we would need an incredibly large computer memory 
capacity and an enormous amount of computer time. In addition the 
computation of just one or two realizations would be of no interest, 
and we would not be able to perform computer experiments with the same 
initial and boundary conditions. Indeed, only statistical 
quantities are of practical meaning in order to describe the 
randomness within the flow. In other words, we first need mean 
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values, then variance and correlations, and also probability density 
functions . 

The method of approach due to Reynolds just described, has been 
applied for cases of non-reacting turbulent flows as early as 1890. 
Mean values for the components of the fluid velocity , the pressure or 
the concentration of diffusing species, and also the variance and 
correlations of the velocity fluctuations and other flow parameters 
have been defined mathematically as well as experimentally (Tennekes 
and Lumley, 1972) . 

Due to the variations and fluctuations of the fluid density which 
are more likely to occur in reacting flows with large temperature 
gradients, it is now common in turbulent combustion studies to change 
the classical definition a little. It appears that the use of mean 
values, weighted by the density p, is more appropriate to handle fluid 
flows with a nature of randomness, irregularity, diffusivity and 
dissipation. 

Let us consider a reacting flow with a velocity component in a 
direction, u a ; density, P ; pressure, P; temperature, T; enthalpy, h; 
viscous stress tensor along a and B directions, T a g; external force, 
such as gravitational acceleration in a direction, g a ; heat flux in a 
direction, q a ; and mass fraction of species s, Y 9 . The four governing 
equations of continuity, momentum, energy, and species can be 
expressed as follows: 


W 
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Continuity Equation: 


dp 9 

— + (pu tt ) =0 a = 1 , 2 , 3 

dt dx a 

Momentum Equation: 

d d 3P 3 t aB 

— (pu a ) + (pu a u 6 ) = pg a - — + 

dt dx 6 ax a dx 6 

Energy Equation: 

d T f PI 9 T f P 

— p h - — + — ugUg + p u 5 h - — 

dt [ \ p 2 / J dx 6 [ l P 

1)1 3 a< l 6 

+ — UgUg = pU^gg - ( w a6 u a) ” 

2 J J 9 X 0 dx fl 

Species Equation: 

9 9 9 

— ( pY 8 ) + (pY 8 u 6 ) (pY 8 u 8 0 ) + W s 

dt dx 6 9 X 0 


(2-1) 


(2-2) 


(2-3) 


(2-4) 


where u 8 0 is the diffusion velocity of species s in & direction. 

In these expressions, viscous stress tensor, T aB ; stress tensor, 
n aB ; and heat flux, q a , can be shown in the following forms: 


9u a 9u 6 ) 2 f 9u k ' 

+ - — U 3, 

9x s 9x tt J 3 ldx k j 


l dx 6 9x tt 
^aS “ ^aB “ T a& 


(2-5) 


(2-6) 


-k t h + P E(¥ s h 8 u 8a ) + q fl 


(2-7) 


where q^ represents the sum of the Dufour heat flux and radiative 
heat flux; k th , the thermal conductivity; and W 8 , the chemical 
reaction rate of species s. 
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To obtain mean values of the thermochemical balance equation, 
density weighted values (or Favre mean values) are used for velocity 
components (u a ) , temperature (T) , mass fraction (Y s ), and enthalpy 
(h) , while the conventional mean values (or Reynolds mean values) are 
kept for pressure (P) and density (/>). 

The conventional mean values are defined as 

A - A + A" (2-8) 

where 7T is the conventional mean value; and () ', the corresponding 
value for fluctuation. The density weighted mean values are defined 
as 

A - A + A' ' (2-8) 

where () is the density weighted mean value; and () ' ' r the 
corresponding value of fluctuation. Detailed mathematical 
expressions for both conventional and density weight mean values are 
discussed in Appendix A. 

The mean value balance equations for continuity, momentum, 


energy, and species concentration can be shown as follows: 
Continuity Equation: 


9 P 9 

— + (pu a ) * 0 

3t 3x a 

Momentum Equation: 


— (/>u a ) + (pu a u B ) 

at ax 6 



ap a 

PQa ~ + ( T a5 

3x a 3x fi 


(2-9) 


( 2 - 10 ) 
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Energy Equation: 
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3P 3 _ ^ 
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3t 9x a 
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— JJ I ii . ~ 

+ T aB u B + T aB u B ) " — 

9x a 

Species Equation: 

a a/ a V 

— ( pY s ) + (pY 8 u a + pY s u a ) 

at ax a 

3 'N-' — ' il 

= (pY s u Stt + pY 8 u' a ) + W s 

3x„ 


( 2 - 11 ) 


( 2 - 12 ) 


where u.„ is the diffusion velocity of species s in a direction. 
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In these manipulations, the following relations have been used: 


T aB * T aB + T afi' 
3a = 3a + 3a / 


(2-13) 

(2-14) 



( 9u a 

a^Bl 


rauj 

9u a 1 

T aB “ u 

+ ‘ 

+ a 

+ " 

laxfi 

ax tt J 


l3x 6 

ax a I 





3ub ) 

2 r 

+ 



+ 

— 



Ux B 

9x a J 

3 k 


9 u k 1 



+ 

u'- 

— 

8 aB 



3x k J 



U = U 

+ it' 




// 


_ 8u k _ 3u k 

a + u 

ax k ax k 


(2-15) 

(2-16) 
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2 
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i 

Ux s 

8x a > 

(X 

3 

8x k J 


_ / 




5 aB 


k t h - k t h + k t h 

( 3T 


(2-17) 

(2-18) 


<Ja 


- Ik 


t h 


+ k 


3x, 


t h 


8T' 

8x, 


+ k 


t h 


8T 

8x, 


/ / \ 


/V 


+ P IY 8 h 8 u 8a 

S 


-D 


+ Y 8 h w # u '' a + Ylh';u sa + Ox + Y'Xu's'a + <** 

8T 

q tt ' = - k th + pi (h 8 Y 8 u 8a + Y 8 h 8 u 8a + u 8tt Y 8 h 8 ) 


(2-19) 


8x 


D 

+ Sa 


( 2 - 20 ) 

For the purpose of numerically solving turbulent reacting flows, 
one has to solve the mean value balance equations of continuity, 
momentum, energy, and species concentration, shown in Equations (2-9) 
to (2-12) , in conjunction with the mean value equation of state 
P = 7R~ <2- 21 > 

Careful examination of Equations (2-9) to (2-12) indicate that 
turbulent diffusion fluxes contain u^[ ujj, u^ Y 8 , and u a h terms which 
are usually much larger than the mean molecular length and time scale 
in the turbulent regime, and also involve correlations between 
fluctuations, and certainly not the mean values only (see Appendix for 
description) • Unless turbulent diffusion fluxes are well defined, 
it is very difficult to obtain a closed form solution numerically. In 
other words, the study of turbulence becomes a closure problem in 
which the modeling of turbulent diffusion fluxes must be pre- 
determined before one can solve the turbulent reacting flow problem 
numerically. 


10 


Ill . Time Scale and Length Scale a£ Tuybujencg 


The mean balance equations of continuity, momentum, energy and 
species, shown in Equations (2-9), (2-10), (2-11) and (2-12), 

respectively, can not be used directly for solving turbulent reacting 
flow problems. These equations involve new quantities the 
turbulent Reynolds stress/diffusion fluxes, such as u^ug, u a Y* and 
u^h?t which are usually much greater than the mean molecular scale in 


the turbulent regime, and also correlations between fluctuations 

which are not mean values only. 

The modeling of the turbulent diffusion fluxes generally follow 
Boussinesq's expression resulting in the following closure 
assumption: 


— •' it 

P u a u B = 


— il v a 
P u a Y 3 


- « Jl 

P u a h 


= Ui 


3u 


__a + &k) 

3x 6 3x a 
/v 


2 

5 aS 

3 


pk + u. x 


yv 
L k 


3u ' 


3x 




ttt fY. 

Sc t 3x 0 
u* 3h 


Pr t 3x tt 


(3-1) 


(3-2) 


(3-3) 


where k - 1/2 u£u£ is the turbulent kinetic energy; u t , the eddy 
viscosity; Sc t , the turbulent Schmidt number; and Pr t , the turbulent 
Prandtl number. It is noted that u t is no£ a property of the fluid. 
In fact, experiments show in particular that Sc t and Pr t are 
independent of the nature of the species, and therefore u t has to be a 
property of the turbulent flow. In other words, u t does depend on a 
length scale and a time scale of the turbulence. This can be shown as 
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(3-4) 


follows with modification from Equation (3-1) : 


Ut 



-<«£«*> l/ * 

p 

r 3u^ 8^ 


A/ /v 

8u a 8 u b ' 



i.8xg 9 x a ; 

m 

k 8x s 8x a ; 

ra 


/ }i 1 / 2 
(U a u 6 ) 


(Length Scale) (Velocity Scale) 

subscript m denotes that the quantity is evaluated where u"u B is a 
maximum value. 

Kolmogorof f suggests that at sufficiently high Reynolds numbers 
there is a range of high wave numbers where the turbulence is 
statistically in equilibrium and uniquely determined by the 
parameters of viscous dissipation of turbulent motions, e, and 
kinematic viscosity, x(= u/p ) (Hinze, 1975). This state of 
equilibrium is universal . This equilibrium range is termed 
"universal" because the turbulence in this range is independent of 
external conditions, and any change in the effective length scale and 
time scale of this turbulence can only be a result of the effect of the 
parameters e and v. The Kolmogorof f length scale, n k , time scale, T k , 

and velocity scale, v k , are represented by: 
fy 3 \ 1/4 




U ) 

f y \ 1 / 2 

i € i 


(3-5) 


(3-6) 


v k =(*0 1/4 (3_7) 

The eddy wave number, k d , where the viscous effects become very strong 
will be of the same order as l/ifc* is defined as: 
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f € \ 1/4 


— - 
’l k 


u 3 y 


(3-8) 


In other words, the viscous effect becomes much stronger than the 
turbulent dissipation when the eddy wave number is much greater than 


Viscous dissipation of turbulent motion is expressed as: 

^ ~ d < ■s^sf 


€ = V 


3x 6 


(3-9) 


8x 6 8x c 



For the case of homogeneous turbulence in which all spatial 
derivatives of mean turbulent quantities become zero, the viscous 
dissipation of turbulent motion becomes: 

(3-10) 

Definitions for the turbulent velocity scale V t , integral length 
scale X j , and Taylor micro-length scale \ t applied to subsequent terms 
yield: 

~~ (3-11) 

(3-12) 


(3-13) 


u a u 6 ^ V t 





, k II it . * 

(u tt u a u B )/v 

9x s Xi 


(3-14) 


where Equation (3-11) shows density-weighted Reynolds stress; 
Equation (3-12) , turbulence production by density-weighted Reynolds 
stress; Equation (3-13) , viscous dissipation of turbulent energy by 
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(3-15) 


Comparison of length scales indicates: 

» x T » i) k viscous dissipation 
where Xj, Xy, and denote integral length, Taylor micro— length, and 
Kolmogoroff length scales, respectively. Equation (3-15) shows a 
cascade transport of turbulent energy in which large inhomogeneous 
eddies transport energy to small homogeneous eddies, which finally 
disappear through viscous dissipation. 
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IV. Examination of Various Turbulence Models 


In the previous section, we have indicated that the diffusion 
fluxes, "3[u and*u£I?obtained from density-weighted mean values 
of the balance equations, require a specification for the closure 
problem. By using Boussinesq's eddy viscosity model, it is evident 
that eddy viscosity K t is not a property of the fluid, but is simply a 
function of a length scale and a time scale of the turbulence. It also 
has been shown in the previous section that there are various length 
scales of turbulence, which come from the inhomogenous large eddy 
integral length scale, Taylor microscale and a homogenous small eddy 
Kolmogorof f length scale, before viscous dissipation becomes 
effective in energy dissipation. 

There are several ways to model the eddy viscosity before one 
can solve the turbulent reacting flow by means of satisfying the 
closure problem. 

(A) Zero Equation Models 

Zero equation models are mostly based on the eddy viscosity 
concept. This concept uses the mathematical convenience of 
retaining the same form of differential equations for laminar and 
turbulent flows and allowing the use of the same solution procedure. 

The first turbulence model proposed, Prandtl's mixing length 
hypothesis, is still widely used. It employs the eddy viscosity 
concept which relates the turbulent transport terms to the local 
gradient of mean flow quantities, namely; 
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// 

u a u B 


- V 


(du 

dx 


a 

+ 


& 


dxj 


(4-1) 


where y t is the eddy viscosity. 

The Prandtl mixing length hypothesis calculates the 
distribution of eddy viscosity by relating it to the mean velocity 
gradient 

A/ 

au sl 


C„L 




fau„ 


ax 6 




(4-2) 


This relation involves a single unknown parameter, the mixing length 
L m whose distribution over the flow field has to be prescribed with an 
empirical information. C tt is a constant. The mixing length model 
has been used for thin shear layers and wall boundary layers 
(Spalding, 1982) . The main drawback with this model is the 
evaluation of L ro for different types of turbulent flows. The 
evaluation of L m becomes difficult for recirculating- flows, three 
dimensional flows, reacting flows, etc. In the already empirical 
specification of the mixing length , it is difficult to incorporate in 
any useful manner , the effects of curvature, buoyancy, rotation, etc. 
The transport and memory effects of turbulence are not accounted for 
in the mixing length model because transport equations for turbulent 


quantities are not involved in this model. 

(B) One Equation Models 

One of the simplest ways to characterize the turbulence is to 
consider the eddy viscosity in terms of two quantities, \ 4 its 
integral length scale and k its kinetic energy (the time scale is 
nothing but \i/k 1/2 ). Eddy viscosity by dimensional analysis 
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becomes : 


(4-3) 


v t - c u Xi k 1/2 

The integral length scale is specified algebraically and hence is 
turbulent flow dependent. The turbulent kinetic energy is derived 


from the Navier-Stokes equations in the one equation model by solving 


its transport properties (Launder and Spalding, 1972) . This 
approach performs only marginally better than the zero equation model 
in terras of the modeling of recirculating flows, three dimensional 
flows, reacting flows, etc. 

(C) Two Equation Models 

In attempts to eliminate the need for specifying the turbulent 
length scale as a function of position throughout the flow, the eddy 
viscosity, shown in Equations (4-2) and (4-3) , can be replaced by the 
following equation: 


~ 2 



In this expression, C u is a scalar constant for isotropic turbulence 
and becomes a vector quantity (no longer a constant value) for the 
marginal extension to cover the cases of anisotropic turbulence, such 
as three-dimensional rotation flows. The eddy viscosity can be 
determined if one can solve two differential equations, one for the 
density weighted mean values of turbulent kinetic energy k and the 
other one for the density weighted mean value of turbulent dissipation 


( . 

The modeled equations for k and e are given below: 
Kinetic Energy (k) Equation: 
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Kinetic Energy Dissipation (c) Equation: 
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u t 3 e 


s t ax tt 


" C £ 2 P ~ 
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(4-6) 




ax 6 


_ ^ au a 

P k - - /> u a u 6 


(4-7) 


(4-8) 


(4-9) 


ax 6 


The coefficients , such as C (1 , C« 2 ' 1 1 , are constants in the sense that 
they are not changed in the modeling calculation. However, these 
constants need to be changed in order to accommodate the effects such 
as curvature, low Reynolds number, near wall effects, etc. Pr t is the 
turbulent Prandtl number which is not species dependent, and is also 


considered a constant in the modeling calculation. 

(D) Two Equation Models Modified by Turbulent Combustion 
The k-f model, shown in the previous section, is modified through 
the interaction of turbulent flow and turbulent combustion. Bray and 
Libby (1976), Bilger (1976), and Jones (1980), and Bray et al., (1981) 
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have devoted particular attention to this problem. The modified 
balance equations for turbulent kinetic energy and dissipation rate 
of turbulence can be shown as follows: 

Modified Turbulent Kinetic Energy Equation: 
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a a 9 < 

— pk + (p u a k) = 
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ax a 3x tt 
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(4-11) 


(4-12) 


(4-13) 


»a - “"'•a 

and Sc t k is the turbulent Schmidt number which is not species 

dependent and is chosen to be 0.9 in Equation (4-12). 

Modified Dissipation Rate Equation of Turbulence: 
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2. T ad is the adiabatic flame temperature and T 0 is the mixture 
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temperature of fuel and oxidizer. 


(E) Reynolds Stress/Flux Models 

The k-e model is based on the eddy viscosity/diffusivity concept 
in which it is assumed that the eddy viscosity/diffusivity is 
isotropic. This means that various Reynolds stress u'^ug and 
diffusion fluxes u^j>"(*" is either Y," or h") are all taking similar 
forms. This assumption is certainly not true for complex flows 
(Jones, 1980) . To overcome this problem, transport equations for the 
various Reynolds stresses and scalar diffusion fluxes must be derived 
from the Navier-Stokes equations. These equations contain higher 
order correlations which have to be approximated by lower order 
correlations to obtain a deterministic set of equations. The 
application of a Reynolds stress/ flux turbulence closure model is 
limited because, especially in flows with a large number of species, 
the n umb er of equations to be solved is large. 

Reynolds Stress Closure Models include a transport equation for 
Reynolds stresses, a transport equation for diffusion scalar fluxes, 
a kinetic energy equation and a kinetic energy dissipation equation. 
These equations are shown as follows (Kollmann and Vandromme, 1979; 
Jones, 1980; Borghi, 1988); 

Transport Equation of Reynolds Stresses: 
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C 2 and C 2 are constant. 


Transport Equation for Diffusion Scalar Fluxes: 
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C # , and C^ 2 are constants. 

Modeling Equation for Fluctuating Velocity Components: 
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where C 8 , C lc , C 2c , C 3c , and C 4e are constants. 
Turbulent Kinetic Energy Equation: 
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where P k is shown in Equation (4-18) , and C, is a constant. 
Turbulent Kinetic Energy Dissipation Equation: 
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(F) Algebraic Stress/Flux Models 

The Reynolds stress/ flux model discussed in Section (4-E) is one 
of the most general turbulence closure models that can be applied to 
almost all kinds of turbulent flows. Disadvantages of this method, 
however, are the large number of differential equations that have to 
be solved and the complexity of the model. To overcome these 


problems. Algebraic Stress Models have been developed in which the 
transport equations for Reynolds stresses and diffusion scalar fluxes 
are simplified to algebraic equations, whereby the model still 
retains most of its basic features, such as for instance the 
anisotropic eddy viscosity/dif fusivity concept (Rodi, 1976). 
Transport Equation for Reynolds Stresses in Algebraic Expression: 
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where Cj and C 2 are the constants. 

Transport Equation for Diffusion Scalar Fluxes in Algebraic 
Expression: 
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where C^ 2 and C^ 2 are constants. 

Modeling Equation for Fluctuating Velocity Components: 
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where C lc , C 2c , C 3c and C 4c are constants. 

(G) Multiple-Scale Turbulence Models 

Turbulence consists of fluctuation motions with a continuous 
spectrum of length and time scales. The largest length scale of 
eddies have the dimension of the flow field, while the smallest length 
scales are related to the diffusive action of the molecular viscosity 
(Tennekes and Lumley, 1972) . Turbulent energy production is related 
to the largest eddies, while dissipation of the turbulent energy is 
primarily related to the smallest eddies. Turbulence closure models 
which employ just a single length scale, therefore, are very 
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simplistic because they neglect the fact that production and 
dissipation of turbulent energy occur at different length scales. 
The fact that single-scale models yield successful predictions of 
many flow fields reflects more that these flows are close to spectral 
equilibrium. In other words, the production of turbulent energy of 

A/ A/ 

the large eddies ( tp) equals the dissipation of the small eddies ( € 7 ) , 
namely 7 ^ = 7 ^ = T, rather than that the single-scale model is adequate 
(Hanjalic et al. r 1980). In view of turbulent energy production 
occurring at large length scales and dissipation of turbulent energy 
mainly taking place at small length scales, there must be a transfer 
region of length scales in which turbulent energy is transferred from 
the largest eddies to the smallest. This transfer can introduce a lag 
phenomenon, namely turbulent energy production and turbulent energy 
dissipation do not both decrease or increase in the same regions of the 
flow (see Figure 1 for details) , as for instance the k-e model implies 
(Hanjalic, et al. , 1980) . Figure 1 shows the energy spectral density 
of turbulence in a shear flow at high Reynolds number as a function of 
the wave number K (1/K is proportional to the length scale) . 

To introduce a model which takes into account the different 
processes at different length scales (or wave number) , Hanjalic et al . 
(1980) divided the energy spectral density distribution into three 
parts, as shown in Figure 1. They assume that production of turbulent 
energy takes place at wave numbers below , characterized by kinetic 
energy^ and energy dissipation 7^, that transfer of turbulent energy 
takes place in the wave number range from Kj to K 2 , characterized by 
kinetic energy k T and energy dissipation < T and that at wave numbers 
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division of the energy spectral density distribution into more than 
three parts is the key feature of the multi-scale-model . The energy 
spectral density range may be divided into as many parts as possible. 
In practice however, a division into three regions appears to be 
sufficient (Fabris et al., 1981). This requires two sets of 
transport equations. The k-e equations for the production region (K 
< K x ) can be written as follows: 

Kinetic Energy Equation for Production Region (K < K x ): 
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Energy Dissipation Equation for Production Region (K < K x ): 
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where P k is shown in Equation (4-18) and C pl and C p2 are constants. 

The k -6 equation for the transfer region (K 2 > K > K x ) can be shown 
as follows: 

Kinetic Energy Equation for the Transfer Region (K 2 > K > K x ) : 
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y Model-Free simulations of Turbulent Reacting Fl9 w -S 


In some cases, the approach based on turbulence models has been 
very successful, and the results obtained by such models have shown 
encouraging agreement in comparison with experimental data. In 
complex chemically reacting flow, however, modeling is extremely 
difficult because of the lack of knowledge of the intricate flow 
dynamics. Also, a turbulence closure optimized for a particular type 
of flow may not be adequate for predicting the flow behavior in other 
configurations. Moreover, since most of the interesting dynamical 
behavior of a flow is modeled a priori, the outcome of numerical 
computations based on these turbulence models cannot substantially 
advance our understanding of turbulent reacting flows (Givi, 1989) . 

Recently, improved efficiency of numerics, storage capability, 
and computational speed have made it possible to solve appropriate 
transport equations of turbulent flows directly, without the need for 
modeling over some limited parameter range (Givi, 1989) . Such 
simulations (defined as Model-Free) , in comparison with calculations 
utilizing turbulence models, have the advantage that the physics of 
the problem is not modeled, but is recovered directly from the 
computed results . The results can be further used to understand many 
important mechanisms of turbulent transport and their direct coupling 
with chemical reactions. 

Model-free simulations consist of solving the time development 
of the detailed, unsteady structures in a turbulent flow field. The 
nonlinear transport equations are solved by means of very accurate 
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numerical methods so that no averaging (Reynolds or Favre average) is 
necessary. Therefore, no turbulence modeling is required. In these 
simulations, data acquisition and statistical sampling of the 
ensemble data are performed by a procedure analogous to laboratory 
experiments. In this respect, the model-free simulations offer the 
advantage that the thermo-fluid parameters of the flow field can be 
easily varied, and the conditions of the numerical simulations are 
more controllable than those in the laboratory experiments (Givi, 
1989) . 

The majority of the approaches in model-free simulations can be 
classified into two general sub-categories: (1) Direct Numerical 
Simulations (DNS), and (2) Large Eddy Simulations (LES) . 

(A) Direct Numerical Simulations 

The main outcome of direct numerical simulations is the 
emergence of various numerical methods that are able to handle large 
variations of length and time scales within a turbulent flow field. 
The n um erical schemes currently in use for direct numerical 
simulations of reacting turbulent flows can be classified into three 
categories: (1) Spectral and pseudospectral methods (Canuto etal., 
1987; Anderson et al., 1984), (2) Lagrangian scheme (Oran and Boris, 
1987) and (3) Finite difference methods (Davis and Moore, 1985) . 

Despite the capability of present day supercomputers in allowing 
the calculations with more than one million grid points, the range of 
length and time scales that can be resolved by direct numerical 
simulations is substantially smaller than those of turbulent flows of 
practical interest. This limits the applications of direct 
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numerical simulations to basic research problems in which the scales 
of the excited modes remain within the band of computationally 
resolved grids (either in physical or spectral domain) . In practice 
this implies that for an accurate simulation, the magnitude of 
viscosity and diffusivity must be large enough to damp out the 
unresolved scales, and the magnitude of the computational time step 
must be kept small enough to capture the correct temporal evolution of 
the flow. These restrictions impose severe limitations for 
practical engineering applications. 

(B) Large Eddy Simulations 

The methodology of large eddy simulations involves the pre- 
filtering of the transport equation by decomposing the transport 
variables into large-scale and sub— grid— scale components. The 
former component is related to the large eddies in the turbulent 
field, whereas the latter is the component containing the small- 
scale fluctuations (Schumann and Friedrich, 1986; 1987) . 

The pre-filtering of the variable j> (x a , t) is performed by means 
of the convolution integral 


<t (Xj, t)> 
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(5-1) 


where F t is an appropriate filter function with a characteristic 
length along a-direction, A a ; <> represents the filtered variable; 
and the integration is over the entire flow field. The remaining 
portion of from the filtered quantity is defined as the sub-grid- 
scale field, and is represented by 
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*' (X a , t) = * (X a , t) - <* (x a , t)> (5-2) 

The Large Eddy Simulation is a combination of the direct 
numerical simulation for the filtered portion of the transport 
variable <*> , and the sub-grid-scale modeling of the small-scale 
component . At first glance, the methodology seems similar to the 
familiar Reynolds averaging. A closer comparison between the two 
averaging procedures, however, indicates the superiority of the pre- 
filtering method in that only the contribution of the small scale 
structures need to be modeled, whereas in Reynolds type averaging 
closures are needed for all the length scales of motion (Givi, 1989) . 
The fact that the small scales of turbulence exhibit a more universal 
character is the main reason to believe that attempting to provide a 
subgrid scale model would be more promising than the previously 
followed procedures based on Reynolds averaging closures. 

An obvious extension of this approach would involve solving a 
transport equation for density-weighted probability functions of the 
subgrid scalars rather than assuming their form. This approach, like 
its counterpart in turbulence modeling, has the advantage that the 
effect of the chemical reations (scalar-scalar correlations) will 
appear in the closed— form. However, models are needed for molecular 
diffusion within the subgrid. 

The rate of progress in this area of direct numerical simulations 
depends on the advance in the fields of (a) numerical methods, (b) 
supercomputer technology, and (c) the mathematical description of 
reactive flow transport. Forseeable developments of advanced 
computational facilities, however, will not be sufficient to relax 
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the restriction of Direct Numerical Simulations to flows having small 
to moderate variations of the characteristic length and time scales 
(Givi, 1989; 1990) . Hence, the boundaries of applicability of Direct 
Numerical Simulations are, and will continue to be, significantly 
restricted. 
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VI Turbulence Model Recommen dation 


As we have already mentioned, the diffusion fluxes require 
specification for the closure problem. It has been shown by the 
Boussinesq model, that the eddy viscosity, v t , is not a property of the 
fluid, but instead a function of a length and time scale. Also the 
length scales evolve from inhomogeneous large eddy integral length 
scales, through Taylor microscale, to homogeneous small eddy 
Kolmogoroff length scale before viscous dissipation effects 
turbulent energy dissipation. 

In the turbulence modeling of the diffusion fluxes, such as u^ug, 
u^Y| and u^h'', through Boussinesq' s eddy viscosity model, there are 
zero equation, one equation, two equation, modified two equation 
models, etc. These models are based on the assumption of isotropic 
turbulence in which the viscous dissipation of turbulent motion is 
simply given by the following expression: 



(6-1) 


In other words, various Reynolds stresses and diffusion fluxes of 
turbulence take similar forms. This assumption is certainly not true 
for complex flows (Jones, 1980) . To overcome this problem, transport 
equations for various Reynolds stresses and scalar diffusion fluxes 
must be derived from the Navier-Stokes equations. The Reynolds 
stress closure models include a transport equation for Reynolds 
stresses, a transport equation for diffusion scalar fluxes, a kinetic 
energy equation, and a kinetic energy dissipation equation. The 
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Reynolds stress/ flux model is one of the most general turbulence 
closure models that can be applied to almost all kinds of turbulent 
flows. Disadvantages of this method, however, are the large number 
of differential equations that have to be solved and the complexity of 
the model. To overcome these problems, algebraic stress models have 
been developed in which the transport equations for Reynolds stresses 
and diffusion scalar fluxes are simplified to algebraic equations, 
whereby the model still retains most of its basic features, such as for 
instance the anisotropic eddy viscosity/dif fusivity concept. 
However, this model has been hardly tested in combustion problems. 

Turbulence consists of fluctuation motions with a continuous 
spectrum of length and time scales. The largest length scale of 
eddies have the dimension of the flow field, while the smallest length 
scales are related to the diffusive action of the molecular viscosity 
(Tennekes and Lumley, 1972) . Turbulent energy production is related 
to the largest eddies, while dissipation of the turbulent energy is 
primarily related to the smallest eddies. The fact that single-scale 
models, employed by turbulence closure models, yield successful 
predictions of many flow fields, reflects more than that these flows 
are close to spectral equilibrium. In other words, the dissipation 
of turbulent energy of the large eddies equals the dissipation of the 
small eddies rather than that the single-scale model is adequate 
(Hanjalic et al., 1980). Multiple-scale turbulence models were 
introduced to take into account the different processes at different 
length scales which includes a transfer region of length scales in 
which turbulent energy is transferred from the largest eddies to the 
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smallest. However, this model has not been tested to treat 
combustion problems. 

In some cases, the results obtained by turbulence models have 
been very successful in comparison with experimental data. In 
complex chemically reacting flows, however, modeling is extremely 
difficult because of lack of knowledge of the intricate flow dynamics. 
Also, a turbulence closure optimized for a particular type of flow may 
not be adequate for predicting the flow behavior in other 
configurations. To avoid these deficiencies caused by turbulence 
modeling, it is suggested to solve appropriate transport equations of 
turbulent flows directly, without the need for modeling over some 
limited parameter range. Model-free simulations have the advantage 
that the physics of the problem is not modeled, but instead taken from 
the computed results. Therefore, the results can be used to 
understand important mechanisms of turbulent transport and their 
direct coupling with chemical reactions. Direct numerical 
simulations and large eddy simulations are a large part of model— free 
simulations. Direct simulations are able to handle large variations 
of length and time scales in a turbulent flow field by the emergence of 
various numerical methods. Direct numerical simulations are limited 
to basic research problems rather than engineering applications 
because the range of length and time scales that can be resolved by 
direct numerical simulations is much smaller than those of turbulent 
flows 

In this report, our major concerns were to examine how the 
available turbulence modelings or the methods of model-free direct 
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simulation of turbulence can be efficiently used for the practical 
engineering application in liquid rocket thrust chambers. With 
these major considerations in mind, the advantages and disadvantages 
of the available turbulence models are summarized below. 

(A) Zero Equation Models 

These models assume isotropic turbulence, and can be applied to 
two-dimensional simple shear flows with mild pressure gradients, mild 
curvature with no flow seperation and/or rotation effects. They can 
be further applied to three-dimensional turbulent boundary layers 
with small cross flows and very mild pressure gradients without 
curvature and rotational effects. They are bad for turbulent flows 
with large cross flows; flows with curvature, rotation and 
separation; pressure and turbulence driven secondary flows; flows 
with abrupt changes in shear rate; shock induced separated flows, etc. 

(B) One Equation Models 

Similar to zero equation models, these models assume isotropic 
turbulence, and can be applied to two dimensional simple shear flows 
without curvature and rotation effects. For the case of two 
dimensional separated flows in a diffuser, they are slightly better 
than the zero equation models. Disadvantages for these models are 
similar to those of the zero equation models. 

(C) Two Equation Models 

These models have been widely tested with good results in 
combustion problems. However, the basic assumption of these models 
are still based on the isotropic eddy viscosity/dif fusivity concept. 
This concept fails in complex flow problems. 
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(D) Two Equation Models Modified for Turbulent Combustion 
Modification of two equation models can greatly improve the 
quality of combustion modeling through the interaction of turbulent 
flow and turbulent combustion. The following examples are given to 
illustrate the modifications: 

(D-l) Modification of Combustion Rate (W) 

To find a model for the combustion rate W is at the core of the 
turbulent combustion modeling. The basic considerations of this 
modeling include: (1) The combustion rate (W i ) , given by chemical 

kinetics, is a molecular rate, and is a highly nonlinear function of 
several variables; and (2) Due to nonlinearity, W is a probability 
function for the encounter of combustion rates between various 
species and is not the value of Wj with the mean values of these 
variables only. 

For a bi-molecular reaction between species A and B (A + B 4 C) one 


obtains 

W A - W B = - k e Y*Y g 


and then 


W, 




- k c Y A Y B = - k c Y A Y B - k c Y A Y 


where k e is a forward reaction rate constant, 
situations can exist, as shown below: 

If Y A and Y b are fluctuating out of phase 
® < o + Y^ < Y^Yg 
If Y a and Y b are fluctuating in phase 

y* a y b ' > 0 4* 

If Y a and Y b prevent reaction 


(6-2) 

(6-3) 

The following 


(6-4) 


(6-5) 
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0 and Y A f 0, Y b * 0 (6-6) 

The effect of temperature fluctuations could play an important 
role when the forward reaction constant k c is a function of T, as in the 
case of combustion. Let us assume: 


k c (T) - k 0 exp 


V 

T , 


(6-7) 


where k 0 and T A are constants. 

The mean reaction rate can be modeled with reference to an 
isothermal bimolecular reaction. For non— isothermal reactions, by 
using Equation (6-7) and expanding exp(- T*/T) in a power series, we 
have mean values of the reaction rate shown as: 
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where the terms P 
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n 

are polynomials of — 

T 


(D-2) Modification of Turbulent Diffusion Fluxes By Combustion 
Further modification of turbulent diffusion fluxes, shown in 
Equation (2-12) , is needed for consideration of the influence of 
combustion on turbulence through the modification of the turbulent 
reaction rate. Turbulent diffusion fluxes u a Y { modified by 
combustion can be shown as follows (Borghi and Dutoya, 1978; Borghi 
and Escudie, 1984; Launder, 1976): 
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In this equation, the term Y " remains because the fluctuation is a 
density weighted fluctuation (that is pY^ - 0 and Y f { f 0). 

(E) Reynolds Stress/Flux Models 

Boussinesq's eddy viscosity/diffustion fluxes models have been 
adopted in the turbulence closure modeling of zero equation, one 
equation, and two equation (including k-c model) models. According 
to Boussinesq's formulation, isotropic turbulence is a foundation in 
the turbulence modeling, but this is not true for complex flows 
(Jones, 1980; Givi, 1989) . To overcome the unappropriate assumption 
of isotropic turbulence for the description of thrust chamber 
reacting flow problems, transport equat ion s^for the various Reynolds 
stresses u([ug, and diffusion fluxes u a h and u a Y [ , must be derived from 
the mean value balance equations of continuity, momentum, energy and 
species concentrations. 

The advantages of adopting Reynolds stress/ flux models in 
turbulent thrust chamber reacting flow problems are that they are 
applicable to almost all turbulent flow problems. There are some 
disadvantages for adopting this model which include: (1) a very 
large number of partial differential equations to be solved, 
especially for multi-component mixtures; (2) the mathematical 
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formulation of this model is very complex; and (3) model is hardly 
tested for combustion reacting flow problems. 

With present day advances in supercomputer technology, the 
former two disadvantages quoted can be easily solved as the new 
generation hardware components and software techniques are 
developed. As to the last disadvantage quoted for the limited number 
of combustion reacting flow tests, it can be improved through the 
endeavor of joint controlled modeling - experiment efforts in reacting 
flow problems. 

In this regard, we strongly recommend the Reynolds stress/flux 
model as a major candidate for the future study in thrust chamber 
turbulent reacting flow problems. 

(F) Algebraic Stress/Flux Models 

To overcome the disadvantages of the large number of 
differential equations that have to be solved and the complexity of 
the mathematical models associated with Reynolds stress/ flux models, 
algebraic models have been developed to simplify the transport 
equations for Reynolds stresses and diffusion scalar fluxes to 
algebraic equations (Rodi, 1976) . The advantages of this model are 
that it is very general in nature, and simpler in mathematical forms 
than that of the Reynolds stress/flux models, but still retains the 
basic features, such as the non-isotropic eddy viscosity/dif fusivity 
concepts. The disadvantages of this model are that it is still very 
complex mathematically in comparison with other models, such as k-f 
models, and it has been little tested with respect to combustion 
problems. 
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(G) Multiple-Scale Turbulence Models 

Length and time scale motions vary constantly in turbulenc. 
Although production and dissipation of turbulent energy occur at 
different length scales, turbulence closure models employ only a 
single length scale (Hanjalic et al. # 1980). By introducing a 
turbulence model with different length scales, one can handle the 
turbulence energy spectral density distribution transition from 
smaller wave numbers (turbulence energy production) to larger wave 
numbers (turbulence energy dissipation) . 

The multiple-scale turbulence models can either be associated 
with eddy viscosity/eddy diffusion concept formulations or with 
Reynolds stress/flux type manipulation. 

(H) Model-Free Simulations 

Instead of seeking turbulence models based on Favre and Reynolds 
averages for obtaining mean values and correlation mean values of 
fluctuations for thermochemical and flow variables of chemically 
reacting flows, one can use a model-free simulation by solving the 
balance equations directly (Givi, 1989; 1990). These model-free 
simulations include Direct Numerical Simulations and Large Eddy 
Simulations. 

The advances in the fields of model— free simulations are 
strongly dependent upon the mathematical description of reactive flow 
transport in addition to the progress in supercomputer technology and 
numerical methods. In other words, direct numerical simulations 
having small variations of characteristic length and time scales will 
still be restricted even with the development of advanced computer 
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facilities (Givi, 1989; 1990). Within a forseeable future, the 
applicability of direct numerical simulation to the turbulent 
reacting flow problems will be mainly in the academic interests, and 
will be significantly restricted for engineering applications. 
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VII. Conclusions 


Turbulence is rotational and three-dimensional. Turbulence is 
characterized by high levels of fluctuating vorticity and random 
fluctuations in nature. Turbulent flows are always dissipative, 
while the random waves are essentially non-dissipative (Tennekes and 
Lumley, 1972) . Turbulence is a continuum phenomenon, and even the 
smallest scales occurring in a turbulent flow are ordinarily far 
larger than any molecular length scale. Turbulence is not a feature 
of fluids but of fluid flows. 

Thermochemical flow parameter and fluctuations of motions of 
length and time scales vary continuously in turbulence. This means a 
turbulent flow has to be considered as an ensemble of random flow 
fields. Indeed, only statistical quantities are of practical 
meaning in order to describe the randomness within the flow. 

Turbulent reacting flow problems can be attacked either by 
model-free direct simulation approaches, or by the use of Favre and 
Reynolds averages for the ensemble of flow fields. Model-free 
simulations include Direct Numerical Simulations and Large Eddy 
Simulations. Model-free simulations will continue to be 
significantly restricted for engineering applications of turbulent 
reacting flow problems because the advancement of supercomputer 
technology and numerical methods alone can not solve the mathematical 
description of reactive flow transport arising from the nature of 
random fluctuations of turbulent flow fields (Givi, 1989; 1990). 
Table 1 shows the classifications of approach used to handle turbulent 
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flow problems including model— free simulation and turbulence 
modeling. 

The other approach is to connsider the variations and 
fluctuations of fluid density, which are more likely to occur in 
reacting flows with large temperature differences in turbulent 
combustion studies, by using mean values weighted by the density to 
handle fluid flows with a nature of randomness, irregularity, 
diffusivity and dissipation. This type of approach can be further 
divided into two large groups, with and without considering the 
Boussinesq expression of closure assumption. 

Turbulent reacting flow problems can not be solved by using mean 
balance equations of continuity, momentum, energy, and species 
concentration. These equations involve new quantities of the 

, ^ /T if 

turbulent Reynold stress/diffusion fluxes. For example, u a ug, u a h, 
and u'^Yj are usually greater than the mean molecular scale, also 
correlations between fluctuations are not mean values only. 

With adopting Boussinesq' s expression of closure assumption, 
one can introduce the concepts of eddy viscosity/eddy diffusion 
fluxes. Turbulence modelings based on Boussinesq' s expression 
include zero equation, one equation, two equation (this contains k-c 
models with and without modifications of combustion rate chemical 
reaction rate, etc. ) , and multiple-scale models associated with eddy 
viscosity. Most of the turbulence modelings based on Boussinesq' s 
expression of closure assumptions have been widely used in 
engineering applications, and also widely tested with good results in 
combustion problems. However, the fundamentals of Boussinesq' s 


44 





expression are built on the assumption of an isotropic eddy 
viscosity/diffusivity concept. This assumption is certainly not 
true for complex flows (Jones, 1980; Givi, 1989) . In particular, the 
isotropic turbulence assumption fails in complex flow problems for 
the description of thrust chamber reacting flow fields. 

Reynolds stress/flux models adopt the transport equations 
derived from the mean values balance equations of continuity , 
momentum, energy and species without using Boussinesq's expression of 
closure assumptions. This makes Reynolds stress/flux models very 
general and applicable to almost all flow problems. It is true that 
this model requires a large number of partial differential equations 
to be solved, especially for multi-component mixtures, in addition to 
the complexity in its mathematical formulation. With the 
advancement in supercomputer technology and numerical methods in the 
foreseeable future, we strongly recommend that the Reynolds 
stress/ flux model be adopted as a major candidate in the future 
development of computation tools for the study of thrust chamber 
turbulent reacting flow problems. 

A workshop was held on the campus of the University of Alabama in 
Huntsville during the time period of April 15 and 17, 1991, in 
connection with one task of a contract between the University and NASA 
Marshall Space Flight Center. The contractual effort was oriented to 
review various turbulence models, used in existing liquid rocket 
thrust chamber flow simulation programs, and other potential 
techniques. The study resulted in the recommendation of the Reynolds 
Stress Model. 
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Appendix A 


Conventional Mean Values (Reynolds Average) 
and Density Weighted Mean Values (Favre Average) of 
Turbulent Reacting Flows 


Due to the variations and fluctuations of the fluid density which 
are more likely to occur in reacting flows with large temperature 
differences, it is rather common, in turbulent combustion studies, to 
change the classical definition of conventional mean values (Reynolds 
average) to density weighted mean values (Favre average) . It appears 
that the use of mean values weighted by density is better able to 
handle fluid flows with randomness, diffusive and dissipative 
irregularity. 

The conventional mean values are defined as: 


A * A +A' (A-l) 

where "(T is the conventional mean value and ()', the corresponding 
value for fluctuation. () can be further defined as: 


A 


1 

lim — 
t->® 2T 


, T 

Adt 
' -T 


(A-2) 


This formulation of conventional mean values can be written as 


A ( x | , t, v) - A (X 4 , t) + A / (Xj, t, v) (A-3) 


where : 


A (Xi, t) 


/CD 

A ( x j , t, x) f ( y ) d k 


( A-4 ) 


f -CD 

The function f ( x) is the probability function of the property x . This 
function is defined such that f ( x ) d x is the probability that the value 


of A at a fixed point (Xj, t) lies in the range from v to v + dv . 
Because v represents all possible values that A can take at a given 

point (x t , t) , it is clear that: 

00 

f ( k ) d k = 1 ( A ~ 5 ) 

The density weighted mean values are defined as: 

A (X i , t, V ) = A(X i , t) + A"(x lf t, V ) (A-6) 

where is the density-weighted mean value, and ()" the 

corresponding value for fluctuations. () can be further defined as 


a/ 

A(Xi, t) 


00 

A(x if t, K)g»(K)dK 


(A-7) 


' -00 

The function g A (y) is the density weighted probability density 
function of A. It is also clear that 


g A ( k ) d k = 1 (A-8) 

* - 00 

There are some close relations between conventional mean values 
and density-weighted mean values . Some useful relations are shown as 
follows: 
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A = 

A(x t , t) f ( v)dv - 
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f ( K ) d K = 
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A = 

A(x i# t)q^(v)dv 

= A 

g* ( v)d/ 

= A 
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A(x j , t)f ( v ) dy = 
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f ( K ) d K = 
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(A— 12) 
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The equivalence of density -weighted mean values, and 
conventional mean values, shown in the Equation (A-15) , is 
particularly important. To give an example, the continuity equation 

is 


dp d 

— + (pu a ) = 0 

at ax tt 

In density -weighted mean values form the equation becomes 


(A-16) 


a p a 

+ (pu a ) = o 

at ax a 

The equivalent conventional mean values form is 


(A- 17 ) 


— + (p u a + p 7 u a 7 ) - 0 ( A-18 ) 

at ax a 

This means that the Favre-averaging continuity equation is simpler in 
form because it contains no terms involving density fluctuations. 

In general, density weighted mean values apply velocity 
components (u a ) , temperature (T) , mass fraction (Y g ) and enthalpy 
(h) , while the conventional mean values are kept for pressure (P) and 
density ( p) . 
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